
####### Replication script for the analysis in the article:                 ###
####### "The institutionalization of a cleavage: How differential treatment ###
####### affects state behavior in the climate negotiations"                 ###
####### by Paula Castro and Marlene Kammerer                                ###
####### published in International Studies Quarterly                        ###
####### April 2021                                                          ###

#### Questions and comments: paula.castro@zhaw.ch
#### Script run using R version 4.0.3 from 2020-10-10


## clear workspace
rm(list= ls())


###--------------------------------------------------------------------------###
## Setting working directory and getting data                               ####
###--------------------------------------------------------------------------###

setwd("~/Dropbox/Differential treatment/PAPERS/Politicization of the climate/Data and analysis/R/REM_2018/Replication")
# setwd("C:/Users/Marlene/Dropbox/Politicization of the climate/Data and analysis/R/REM_2018/Replication")

load("data_type.RData")
load("data_rate_coop.RData")
load("data_rate_confl.RData")
load("indep_vars_annual.RData")

## Loading necessary packages

library(ggplot2)
library(texreg)
library(dplyr)
library(gridExtra)
library(ggeffects)
library(igraph)


###--------------------------------------------------------------------------###
##   Figures 1 and 2: Networks of cooperative and of conflictual events    #####
###--------------------------------------------------------------------------###

### Subsetting event sequences: only cooperative interactions and only conflictual interactions, and only sender and target

ENBcoop <- subset(data_type, data_type$cooperation == 1, select = c(sender, target))
ENBconfl <- subset(data_type, data_type$cooperation == 0, select = c(sender, target))

### Creating networks/graphs/matrices

ENBcoop_net <- graph.data.frame(ENBcoop, directed=TRUE)
ENBcoop_Edge <- get.edgelist(ENBcoop_net)
ENBcoop_Mat<-get.adjacency(ENBcoop_net, sparse=FALSE)

ENBconfl_net <- graph.data.frame(ENBconfl, directed=TRUE)
ENBconfl_Edge <- get.edgelist(ENBconfl_net)
ENBconfl_Mat<-get.adjacency(ENBconfl_net, sparse=FALSE)

igraph.options(sparsematrices=FALSE)


### Network of cooperative events: dichotomize with threshold = 50 (one tie means at least 50 cooperative interactions)

allcoop_50 <- matrix (0 , nrow = nrow (ENBcoop_Mat) , ncol = ncol (ENBcoop_Mat)) 

for (i in 1:nrow(allcoop_50)) {
  for (j in 1:ncol(allcoop_50)) {
    if (ENBcoop_Mat [i,j] < 50) {
      allcoop_50 [i,j] <- 0 }
    if (ENBcoop_Mat [i,j] >= 50) {
      allcoop_50[i,j] <- 1
    }}}

row.names(allcoop_50) = row.names(ENBcoop_Mat)
colnames(allcoop_50) = colnames(ENBcoop_Mat)

### Create network graph of cooperative negotiation interactions

detach(package:igraph)
library(sna)
library(GGally)
library(network)
library(ggplot2)
library(intergraph)
library(scales)

net_all50 <- network(allcoop_50, directed = TRUE) # create network object
delete.vertices(net_all50, which(degree(net_all50 )==0)) # exclude isolates
plot(net_all50) # create simple plot to check if all is correct
mat.all50 <- as.sociomatrix(net_all50) # create matrix

# Check correlation between full and dichotomized network
qaptest_50 <- qaptest(list(ENBcoop_Mat, allcoop_50), gcor, g1 = 1, g2 = 2)
summary(qaptest_50$test) # correlation coefficient = 0.78

table(data_type$sender, data_type$annexi_sender)

# Add attribute: Group membership
annex <- data_type %>% select(sender, annexi_sender) %>%
  rename(Country = sender, annexi = annexi_sender) %>%
  mutate(Group = if_else(annexi > 0, "Annex I", "non-Annex I")) %>%
  distinct() %>%
  filter(!(Country == "Slovenia" & Group == "non-Annex I")) %>%
  filter(!(Country == "Croatia" & Group == "non-Annex I"))
# Croatia and Slovenia change groups - they become members of Annex I when they join the EU. For the graph, we therefore classify them as "Annex I" countries.

net_all50%v%"Group" <- annex$Group[match(net_all50%v%"vertex.names", annex$Country)]

# Graph 
dev.off()

ggnet2(net_all50, alpha = 0.75, label = TRUE, label.size = 2.2, vjust = -1.2, mode = "kamadakawai",  
       size = "indegree", size.min = 1, size.cut = 5, size.legend = "Indegree", max_size = 4,
       shape = "Group", shape.legend = "Group membership",
       shape.palette = c("non-Annex I" = 16, "Annex I" = 17, "mixed" = 18),
       color = "Group", color.legend = "Group membership",
       color.palette = c("non-Annex I" = "grey15", "Annex I" = "grey65", "mixed" = "grey50"),
       legend.size = 7.5, legend.position = "right") +
  theme(legend.spacing.y = unit(0.2, 'cm'), 
        legend.key.size = unit(1.2, 'lines'),
        legend.title = element_text(face = "bold"))

ggsave("Figure_2.pdf", width = 6.5, height = 4,  device = "pdf", dpi = 300, units = "in")


### Network of conflictual events: dichotomize with threshold = 15 (one tie means at least 15 conflictual interactions)

allconfl_15 <- matrix (0 , nrow = nrow (ENBconfl_Mat) , ncol = ncol (ENBconfl_Mat)) 

for (i in 1:nrow( allconfl_15)) {
  for (j in 1:ncol( allconfl_15)) {
    if (ENBconfl_Mat [i,j] < 15) {
      allconfl_15 [i,j] <- 0 }
    if (ENBconfl_Mat [i,j] >= 15) {
      allconfl_15[i,j] <- 1
    }}}

row.names(allconfl_15) = row.names(ENBconfl_Mat)
colnames(allconfl_15) = colnames(ENBconfl_Mat)

### Create network graph of conflictual negotiation interactions

net_confl15 <- network(allconfl_15, directed = TRUE) # create network object
delete.vertices(net_confl15, which(degree(net_confl15 )==0)) # exclude isolates
plot(net_confl15) # create simple plot to check if all is correct
mat.confl15 <- as.sociomatrix(net_confl15) # create matrix

# Check correlation between full and dichotomized network
qaptest_15 <- qaptest(list(ENBconfl_Mat, allconfl_15), gcor, g1 = 1, g2 = 2)
summary(qaptest_15$test) # correlation coefficient = 0.78

# Attach attribute: group membership
net_confl15%v%"Group" <- annex$Group[match(net_confl15%v%"vertex.names", annex$Country)]

# Graph
dev.off()
ggnet2(net_confl15, alpha = 0.75, label = TRUE, label.size = 2.2, vjust = -1.2, mode = "kamadakawai", 
       size = "indegree", size.min = 1, size.cut = 5, size.legend = "Indegree", max_size = 4,
       shape = "Group", shape.legend = "Group membership",
       shape.palette = c("non-Annex I" = 16, "Annex I" = 17, "mixed" = 18),
       color = "Group", color.legend = "Group membership",
       color.palette = c("non-Annex I" = "grey15", "Annex I" = "grey65", "mixed" = "grey50"),
       legend.size = 7, legend.position = "right") +
  theme(legend.spacing.y = unit(0.2, 'cm'), 
        legend.key.size = unit(1.2, 'lines'),
        legend.title = element_text(face = "bold"))

ggsave("Figure_3.pdf", width = 6.5, height = 4,  device = "pdf", dpi = 300, units = "in")


###--------------------------------------------------------------------------###
##          Type model: Likelihood of cooperation                          #####
###--------------------------------------------------------------------------###

# Create all_mitigation variable
data_type$all_mitigation <- "Other topics"
data_type$all_mitigation[data_type$topic=="Mitigation"] <- "Mitigation"
data_type$all_mitigation[data_type$topic=="Mitigation Annex1"] <- "Mitigation"
data_type$all_mitigation[data_type$topic=="Mitigation non-Annex1"] <- "Mitigation"
data_type$all_mitigation <- as.factor(data_type$all_mitigation)
table(data_type$topic, data_type$all_mitigation)

# Set "Other topics" as reference category
data_type <- within(data_type, all_mitigation <- relevel(all_mitigation, ref = "Other topics"))
levels(data_type$all_mitigation)

names(data_type)

## Table 1: Main results #####

# Model 1 (binomial model without network statistics)
typemodel.1     <-  glm (cooperation ~ 
                           + same_annex_factor*all_mitigation
                         + log_pop_comp_sender + log_pop_comp_target
                         + log_gdpcap_comp_sender + log_gdpcap_comp_target
                         + log_co2cap_comp_sender + log_co2cap_comp_target
                         + democracy_comp_sender + democracy_comp_target
                         + nd_gain_sender + nd_gain_target
                         + forest_area_comp_sender + forest_area_comp_target
                         + fossilrents_comp_sender + fossilrents_comp_target
                         + English_sender + English_target
                         + log_trade
                         + log_aid
                         + same_region
                         + same_coalition_noG77
                         + coalition_member
                         + interv_sender + interv_target
                         + yeardummy
                         , family = "binomial", data = data_type)


# Model 2 (binomial model with network statistics)
typemodel.2     <-  glm (cooperation ~ 
                           + sender_outdegree + target_indegree
                         + reciprocity
                         + triad_fof + triad_foe + triad_eoe + triad_eof
                         + sender_similarity + target_similarity
                         + same_annex_factor*all_mitigation
                         + log_pop_comp_sender + log_pop_comp_target
                         + log_gdpcap_comp_sender + log_gdpcap_comp_target
                         + log_co2cap_comp_sender + log_co2cap_comp_target
                         + democracy_comp_sender + democracy_comp_target
                         + nd_gain_sender + nd_gain_target
                         + forest_area_comp_sender + forest_area_comp_target
                         + fossilrents_comp_sender + fossilrents_comp_target
                         + English_sender + English_target
                         + log_trade
                         + log_aid
                         + same_region
                         + same_coalition_noG77
                         + coalition_member
                         + yeardummy
                         , family = "binomial", data = data_type)


# Model 3 (Linear Probability Model without network statistics)
typemodel.3     <-  lm (cooperation ~ 
                          + same_annex_factor*all_mitigation
                        + log_pop_comp_sender + log_pop_comp_target
                        + log_gdpcap_comp_sender + log_gdpcap_comp_target
                        + log_co2cap_comp_sender + log_co2cap_comp_target
                        + democracy_comp_sender + democracy_comp_target
                        + nd_gain_sender + nd_gain_target
                        + forest_area_comp_sender + forest_area_comp_target
                        + fossilrents_comp_sender + fossilrents_comp_target
                        + English_sender + English_target
                        + log_trade
                        + log_aid
                        + same_region
                        + same_coalition_noG77
                        + coalition_member
                        + interv_sender + interv_target
                        + yeardummy
                        , data = data_type)


# Model 4 (full Linear Probability Model with network statistics)
typemodel.4     <-  lm (cooperation ~ 
                          + sender_outdegree + target_indegree
                        + reciprocity
                        + triad_fof + triad_foe + triad_eoe + triad_eof
                        + sender_similarity + target_similarity
                        + same_annex_factor*all_mitigation
                        + log_pop_comp_sender + log_pop_comp_target
                        + log_gdpcap_comp_sender + log_gdpcap_comp_target
                        + log_co2cap_comp_sender + log_co2cap_comp_target
                        + democracy_comp_sender + democracy_comp_target
                        + nd_gain_sender + nd_gain_target
                        + forest_area_comp_sender + forest_area_comp_target
                        + fossilrents_comp_sender + fossilrents_comp_target
                        + English_sender + English_target
                        + log_trade
                        + log_aid
                        + same_region
                        + same_coalition_noG77
                        + coalition_member
                        + yeardummy
                        , data = data_type)


# Results table:
htmlreg(list(typemodel.1, typemodel.2, typemodel.3, typemodel.4), 
        file = "Table_1.doc", inline.css = FALSE, digits = 3,
        doctype = TRUE, html.tag = TRUE, head.tag = TRUE, body.tag = TRUE, single.row = TRUE)


# Generate predicted probabilities
predict.tm1 <- ggpredict(typemodel.1, terms = c("all_mitigation", "same_annex_factor"))
predict.tm2 <- ggpredict(typemodel.2, terms = c("all_mitigation", "same_annex_factor"))
predict.tm3 <- ggpredict(typemodel.3, terms = c("all_mitigation", "same_annex_factor"))
predict.tm4 <- ggpredict(typemodel.4, terms = c("all_mitigation", "same_annex_factor"))


## Figure 4: Interaction plots of models 1 and 2 #####

# Relevel the categories of all_mitigation, so that the graphs appear in the desired order
predict.tm1 <- within(predict.tm1, x <- relevel(x, ref = "Mitigation"))
predict.tm1 <- predict.tm1[order(predict.tm1$x, predict.tm1$group),]

predict.tm2 <- within(predict.tm2, x <- relevel(x, ref = "Mitigation"))
predict.tm2 <- predict.tm2[order(predict.tm2$x, predict.tm2$group),]

# Individual graphs
gplot.tm1 <- ggplot(predict.tm1, aes(x, predicted, group = group, shape = group)) + geom_line() + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Negotiation topic") + 
  ylab("Predicted probability of cooperation") +
  ggtitle("(a): Without network statistics, logit (Model 1)") + 
  guides(shape = guide_legend(title = "Group membership")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5, margin=margin(0,0,10,0)),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))

gplot.tm2 <- ggplot(predict.tm2, aes(x, predicted, group = group, shape = group)) + geom_line() + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Negotiation topic") + 
  ylab("Predicted probability of cooperation") +
  ggtitle("(b): With network statistics, logit (Model 2)") + 
  guides(shape = guide_legend(title = "Group membership")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5, margin=margin(0,0,10,0)),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Overall figure
plots_fig4 <- list(gplot.tm1, gplot.tm2)
margin = theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))
fig_4 <- grid.arrange(grobs = lapply(plots_fig4, "+", margin))
fig_4

ggsave("Figure_4.ps", fig_4, device = "eps", height = 8, width = 6.5, dpi = 300, units = "in")


## Table A2: Robustness checks for the type model #####

# Create criticaltopic variable
data_type$criticaltopic <- "Other topics"
data_type$criticaltopic[data_type$all_mitigation == "Mitigation"] <- "Critical topic"
data_type$criticaltopic[data_type$topic == "Principles"] <- "Critical topic"
data_type$criticaltopic[data_type$topic == "Content of new agreement"] <- "Critical topic"
data_type$criticaltopic <- as.factor(data_type$criticaltopic)
table(data_type$topic, data_type$criticaltopic)

# Set "Other topics" as reference category
data_type <- within(data_type, criticaltopic <- relevel(criticaltopic, ref = "Other topics"))
levels(data_type$criticaltopic)

# Model A1 (criticaltopic instead of all_mitigation)
typemodel.A1     <-  glm (cooperation ~ 
                           + sender_outdegree + target_indegree
                         + reciprocity
                         + triad_fof + triad_foe + triad_eoe + triad_eof
                         + sender_similarity + target_similarity
                         + same_annex_factor*criticaltopic
                         + log_pop_comp_sender + log_pop_comp_target
                         + log_gdpcap_comp_sender + log_gdpcap_comp_target
                         + log_co2cap_comp_sender + log_co2cap_comp_target
                         + democracy_comp_sender + democracy_comp_target
                         + nd_gain_sender + nd_gain_target
                         + forest_area_comp_sender + forest_area_comp_target
                         + fossilrents_comp_sender + fossilrents_comp_target
                         + English_sender + English_target
                         + log_trade
                         + log_aid
                         + same_region
                         + same_coalition_noG77
                         + coalition_member
                         + yeardummy
                         , family = "binomial", data = data_type)


# Model A2 (total emissions instead of per capita emissions)
typemodel.A2     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2_comp_sender + log_co2_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition_noG77
                          + coalition_member
                          + yeardummy
                          , family = "binomial", data = data_type)


# Model A3 (adding UN voting similarity)
typemodel.A3     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2cap_comp_sender + log_co2cap_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + un_voting
                          + same_coalition_noG77
                          + coalition_member
                          + yeardummy
                          , family = "binomial", data = data_type)


# Model A4 (using same_coalition instead of same_coalition_noG77)
typemodel.A4     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2cap_comp_sender + log_co2cap_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition
                          + coalition_member
                          + yeardummy
                          , family = "binomial", data = data_type)


# Model A5 (using meeting fixed effects instead of year fixed effects)
typemodel.A5     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2cap_comp_sender + log_co2cap_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition_noG77
                          + coalition_member
                          + meeting
                          , family = "binomial", data = data_type)


# Model A6 (excluding coalitions)
data_nocoalitions <- subset(data_type, data_type$coalition == 0)

typemodel.A6     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2cap_comp_sender + log_co2cap_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition_noG77
                          + yeardummy
                          , family = "binomial", data = data_nocoalitions)


# Model A7 (no completed data)
typemodel.A7     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_population_sender + log_population_target
                          + log_gdpcap_sender + log_gdpcap_target
                          + log_co2cap_sender + log_co2cap_target
                          + democracy_sender + democracy_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_sender + forest_area_target
                          + fossil_rents_sender + fossil_rents_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition_noG77
                          + coalition_member
                          + yeardummy
                          , family = "binomial", data = data_type)


# Model A8 (deleting 2 years of initial observations to avoid bias in network statistics)

summary(data_type$eventTime) # in total we have 6864 time periods
data_nostart <- subset(data_type, data_type$eventTime > 730)

typemodel.A8     <-  glm (cooperation ~ 
                            + sender_outdegree + target_indegree
                          + reciprocity
                          + triad_fof + triad_foe + triad_eoe + triad_eof
                          + sender_similarity + target_similarity
                          + same_annex_factor*all_mitigation
                          + log_pop_comp_sender + log_pop_comp_target
                          + log_gdpcap_comp_sender + log_gdpcap_comp_target
                          + log_co2cap_comp_sender + log_co2cap_comp_target
                          + democracy_comp_sender + democracy_comp_target
                          + nd_gain_sender + nd_gain_target
                          + forest_area_comp_sender + forest_area_comp_target
                          + fossilrents_comp_sender + fossilrents_comp_target
                          + English_sender + English_target
                          + log_trade
                          + log_aid
                          + same_region
                          + same_coalition_noG77
                          + coalition_member
                          + yeardummy
                          , family = "binomial", data = data_nostart)


# Print table
htmlreg(list(typemodel.A1, typemodel.A2, typemodel.A3, typemodel.A4, typemodel.A5, typemodel.A6, typemodel.A7, typemodel.A8), 
        file = "Table_A2.doc", inline.css = FALSE, doctype = TRUE, digits = 3,
        html.tag = TRUE, head.tag = TRUE, body.tag = TRUE, single.row = TRUE)


## Figure A9: Interaction plots of further models #####

# Predicted probabilities
predict.tmA1 <- ggpredict(typemodel.A1, terms = c("criticaltopic", "same_annex_factor"))
predict.tmA2 <- ggpredict(typemodel.A2, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA3 <- ggpredict(typemodel.A3, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA4 <- ggpredict(typemodel.A4, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA5 <- ggpredict(typemodel.A5, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA6 <- ggpredict(typemodel.A6, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA7 <- ggpredict(typemodel.A7, terms = c("all_mitigation", "same_annex_factor"))
predict.tmA8 <- ggpredict(typemodel.A8, terms = c("all_mitigation", "same_annex_factor"))

# Relevel
predict.tm3 <- within(predict.tm3, x <- relevel(x, ref = "Mitigation"))
predict.tm4 <- within(predict.tm4, x <- relevel(x, ref = "Mitigation"))
predict.tmA1 <- within(predict.tmA1, x <- relevel(x, ref = "Critical topic"))
predict.tmA2 <- within(predict.tmA2, x <- relevel(x, ref = "Mitigation"))
predict.tmA3 <- within(predict.tmA3, x <- relevel(x, ref = "Mitigation"))
predict.tmA4 <- within(predict.tmA4, x <- relevel(x, ref = "Mitigation"))
predict.tmA5 <- within(predict.tmA5, x <- relevel(x, ref = "Mitigation"))
predict.tmA6 <- within(predict.tmA6, x <- relevel(x, ref = "Mitigation"))
predict.tmA7 <- within(predict.tmA7, x <- relevel(x, ref = "Mitigation"))
predict.tmA8 <- within(predict.tmA8, x <- relevel(x, ref = "Mitigation"))


list_predict_tm3_tm4 <- list(predict.tm3, predict.tm4)
list_predict_A1_A8 <- list(predict.tmA1, predict.tmA2, predict.tmA3, predict.tmA4, predict.tmA5, predict.tmA6, predict.tmA7, predict.tmA8)

# Graphs
for(i in 1:2){
  predict_i <- list_predict_tm3_tm4[[i]]
  predict_i <- predict_i[order(predict_i$x, predict_i$group),]
  gplot.tm_i <- ggplot(predict_i, aes(x, predicted, group = group, shape = group)) + geom_line() + 
    geom_point(size=2) + 
    geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
    xlab("Negotiation topic") + 
    ylab("Predicted probability of cooperation") +
    ggtitle(paste0("Model ", i+2)) + 
    guides(shape = guide_legend(title = "Group membership")) +
    theme_bw() +
    theme(
      plot.title = element_text(size=10, face="bold", hjust = 0.5),
      axis.title.x = element_text(size=8, face="bold"),
      axis.title.y = element_text(size=8, face="bold"),
      legend.text = element_text(size=8),
      legend.title = element_text(size=8, face="bold"))
  assign(paste0("gplot.tm", i+2), gplot.tm_i) 
}

for(i in 1:8){
  predict_i <- list_predict_A1_A8[[i]]
  predict_i <- predict_i[order(predict_i$x, predict_i$group),]
  gplot.tm_i <- ggplot(predict_i, aes(x, predicted, group = group, shape = group)) + geom_line() + 
    geom_point(size=2) + 
    geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
    xlab("Negotiation topic") + 
    ylab("Predicted probability of cooperation") +
    ggtitle(paste0("Model A", i)) + 
    guides(shape = guide_legend(title = "Group membership")) +
    theme_bw() +
    theme(
      plot.title = element_text(size=10, face="bold", hjust = 0.5),
      axis.title.x = element_text(size=8, face="bold"),
      axis.title.y = element_text(size=8, face="bold"),
      legend.text = element_text(size=8),
      legend.title = element_text(size=8, face="bold"))
  assign(paste0("gplot.tmA", i), gplot.tm_i) 
}


# Putting all graphs into two pdfs, including also the plots for models 3 and 4 in the main text of the manuscript:
plots_figA9_a <- list(gplot.tm3, gplot.tm4, gplot.tmA1, gplot.tmA2, gplot.tmA3, gplot.tmA4)
plots_figA9_b <- list(gplot.tmA5, gplot.tmA6, gplot.tmA7, gplot.tmA8)
margin = theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))

dev.off()
fig_A9_a <- grid.arrange(grobs = lapply(plots_figA9_a, "+", margin), nrow = 3, ncol = 2)
fig_A9_b <- grid.arrange(grobs = lapply(plots_figA9_b, "+", margin), nrow = 2, ncol = 2)

ggsave("Figure_A9a.ps", fig_A9_a, device = "eps", height = 12, width = 12, dpi = 300)
ggsave("Figure_A9b.ps", fig_A9_b, device = "eps", height = 8, width = 12, dpi = 300)


## Figure A10: Interaction plots under different control variable values ####

# Using main model (typemodel.2), we predict probabilities of cooperation for different years, 
# assuming that both countries in a pair are members of the same coalition and are located in the same region

predict.tm2.variousyears <- ggpredict(typemodel.2, terms = c("all_mitigation", "same_annex_factor", "yeardummy [1995, 2000, 2002, 2005, 2007, 2009, 2011, 2013]"), condition = c(same_coalition_noG77 = 1, same_region = 1))
predict.tm2.variousyears

# Relevel and reorder
predict.tm2.variousyears <- within(predict.tm2.variousyears, x <- relevel(x, ref = "Mitigation"))
predict.tm2.variousyears <- predict.tm2.variousyears[order(predict.tm2.variousyears$x, predict.tm2.variousyears$group),]

# Graph
gplot.tm2.variousyears <- ggplot(predict.tm2.variousyears, aes(x, predicted, group = group, shape = group)) + geom_line() + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  facet_wrap(~facet) +
  xlab("Negotiation topic") + 
  ylab("Predicted probability of cooperation") +
  ggtitle("Model 2, at various years, assuming pairs are in the same coalition and region") + 
  guides(shape = guide_legend(title = "Group membership")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))
ggsave(filename="Figure_A10.pdf", plot = gplot.tm2.variousyears, width = 20, height = 14, units = "cm")


## Figures A11 - A12: Sensitivity analysis of omitted variable bias ####

library(sensemakr)
# The sensitivity analysis works only with an OLS model that has a binary treatment variable.
# I will therefore separate my original main variable (same_annex_factor) into 2 separate dummies, 
# so that I can run the sensitivity analysis on each of the dummies separately.
# I will also create the interactions manually, so that I can test their sensitivity as well.

data_type$same_annex_v2 <- as.factor(data_type$same_annex)

data_type$both_annexi <- 0
data_type$both_annexi[data_type$same_annex_v2 == 2] <- 1

data_type$both_nonannexi <- 0
data_type$both_nonannexi[data_type$same_annex_v2 == 1] <- 1

class(data_type$both_annexi)
class(data_type$both_nonannexi)

data_type$mitigation <- 1
data_type$mitigation[data_type$all_mitigation == "Other topics"] <- 0
class(data_type$mitigation)

data_type$bothNAI.x.mitigation <- data_type$both_nonannexi * data_type$mitigation
data_type$bothAI.x.mitigation <- data_type$both_annexi * data_type$mitigation

# Renaming some explanatory variables that will be used for the sensitivity analysis:
data_type$same.coalition <- data_type$same_coalition_noG77

# Regression model
typemodel.4.dummies     <-  lm (cooperation ~ 
                                  + sender_outdegree + target_indegree
                                + reciprocity
                                + triad_fof + triad_foe + triad_eoe + triad_eof
                                + sender_similarity + target_similarity
                                + both_annexi
                                + both_nonannexi
                                + bothAI.x.mitigation
                                + bothNAI.x.mitigation
                                + mitigation
                                + log_pop_comp_sender + log_pop_comp_target
                                + log_gdpcap_comp_sender + log_gdpcap_comp_target
                                + log_co2cap_comp_sender + log_co2cap_comp_target
                                + democracy_comp_sender + democracy_comp_target
                                + nd_gain_sender + nd_gain_target
                                + forest_area_comp_sender + forest_area_comp_target
                                + fossilrents_comp_sender + fossilrents_comp_target
                                + English_sender + English_target
                                + log_trade
                                + log_aid
                                + same_region
                                + same.coalition
                                + coalition_member
                                + yeardummy
                                , data = data_type)

screenreg(list(typemodel.2, typemodel.4, typemodel.4.dummies))
# Although it is parameterized differently, the results are exactly the same as in typemodel.4.


# run sensemakr for sensitivity analysis: using two different control variables as benchmark for comparison 
sensitivity.bothAI.samecoalition <- sensemakr(typemodel.4.dummies, treatment = "both_annexi",
                                              benchmark_covariates = "same.coalition1",
                                              kd = c(1, 5, 10))

sensitivity.bothNAI.samecoalition <- sensemakr(typemodel.4.dummies, treatment = "both_nonannexi",
                                               benchmark_covariates = "same.coalition1",
                                               kd = c(1, 5, 10))


sensitivity.bothAI.reciprocity <- sensemakr(typemodel.4.dummies, treatment = "both_annexi",
                                            benchmark_covariates = "reciprocity",
                                            kd = c(1, 5, 10))

sensitivity.bothNAI.reciprocity <- sensemakr(typemodel.4.dummies, treatment = "both_nonannexi",
                                             benchmark_covariates = "reciprocity",
                                             kd = c(1, 5, 10))

# explanation of results
summary(sensitivity.bothAI.samecoalition)
summary(sensitivity.bothNAI.samecoalition)

summary(sensitivity.bothAI.reciprocity)
summary(sensitivity.bothNAI.reciprocity)

# plot omitted variable bias contour of t-value

pdf("Figure_A11.pdf", pointsize = 10, width = 10, height = 5)
par(mfrow=c(1,2), oma=c(2.5,2.5,2.5,2.5), mai=c(2,2,1,1), cex=.7) # set the plotting area into a 2*2 array, increases margins between plots, reduces font size
plot(sensitivity.bothAI.reciprocity, sensitivity.of = "t-value", label.text = F, label.bump.y = 0.01, main = "Contour plot of omitted variable bias, treatment = both Annex I")
text(0.035, 0.05, "1x reciprocity \n(coef=0.15, t=22.01)")
text(0.04, 0.12, "5x reciprocity \n(coef=0.11, t=16.77)")
text(0.06, 0.22, "10x reciprocity \n(coef=0.06, t=9.59)")
plot(sensitivity.bothNAI.reciprocity, sensitivity.of = "t-value", label.text = F, label.bump.y = 0.01, main = "Contour plot of omitted variable bias, treatment = both non-Annex I")
text(0.035, 0.05, "1x reciprocity \n(coef=0.19, t=31.62)")
text(0.04, 0.12, "5x reciprocity \n(coef=0.17, t=29.69)")
text(0.06, 0.22, "10x reciprocity \n(coef=0.15, t=27.10)")
dev.off()


pdf("Figure_A12.pdf", pointsize = 10, width = 10, height = 5)
par(mfrow=c(1,2), oma=c(2.5,2.5,2.5,2.5), mai=c(2,2,1,1), cex=.7) # set the plotting area into a 2*2 array, increases margins between plots, reduces font size
plot(sensitivity.bothAI.samecoalition, sensitivity.of = "t-value", label.text = F, label.bump.x = 0.01, main = "Contour plot of omitted variable bias, treatment = both Annex I")
text(0.06, 0.055, "1x same.coalition \n(coef=0.14, t=20.44)")
text(0.21, 0.06, "5x same.coalition \n(coef=0.07, t=8.88)")
text(0.42, 0.065, "10x same.coalition \n(coef=-0.06, t=-6.41)")
plot(sensitivity.bothNAI.samecoalition, sensitivity.of = "t-value", label.text = F, label.bump.x = 0.01, main = "Contour plot of omitted variable bias, treatment = both non-Annex I")
text(0.12, 0.052, "1x, 5x, 10x same.coalition \n(coef=0.19, t=32; coef=0.18, t=30; coef=0.17, t=27)")
dev.off() 



###--------------------------------------------------------------------------###
##          Rate (survival) model: Frequency of cooperation                #####
###--------------------------------------------------------------------------###

# Model 5 (cooperative events, year dummies; no network variables)
ratemodel.5 <- glm(eventDummy ~ 
                     same_annex_factor*yeardummy +
                     log_pop_comp_sender + log_pop_comp_target +
                     log_gdpcap_comp_sender + log_gdpcap_comp_target +
                     log_co2cap_comp_sender + log_co2cap_comp_target + 
                     democracy_comp_sender + democracy_comp_target +
                     nd_gain_sender + nd_gain_target +
                     forest_area_comp_sender + forest_area_comp_target +
                     fossilrents_comp_sender + fossilrents_comp_target +
                     English_sender + English_target +
                     log_trade +
                     log_aid +
                     same_region +
                     same_coalition_noG77 +
                     coalition_member +
                     interv_sender + interv_target +
                     yeardummy,
                   family = binomial(link = logit), data = data_rate_coop)

# Predicted probabilities: 
predict.rm5 <- ggpredict(ratemodel.5, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rm5 <- ggplot(predict.rm5, aes(x, predicted, group = group, shape = group, linetype = group)) + 
  geom_line(size = .6) + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a cooperative event happens") +
  ggtitle("Rate Model 5") + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model 6 (cooperative events, year dummies, all variables)
ratemodel.6 <- glm(eventDummy ~ 
                     sender_outdegree + target_indegree +
                     reciprocity + 
                     triad_fof + triad_foe + triad_eoe + triad_eof +
                     sender_similarity + target_similarity +
                     same_annex_factor*yeardummy +
                     log_pop_comp_sender + log_pop_comp_target +
                     log_gdpcap_comp_sender + log_gdpcap_comp_target +
                     log_co2cap_comp_sender + log_co2cap_comp_target +
                     democracy_comp_sender + democracy_comp_target +
                     nd_gain_sender + nd_gain_target +
                     forest_area_comp_sender + forest_area_comp_target +
                     fossilrents_comp_sender + fossilrents_comp_target +
                     English_sender + English_target +
                     log_trade +
                     log_aid +
                     same_region +
                     same_coalition_noG77 +
                     coalition_member +
                     yeardummy,
                   family = binomial(link = logit), data = data_rate_coop)

# Predicted probabilities: 
predict.rm6 <- ggpredict(ratemodel.6, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rm6 <- ggplot(predict.rm6, aes(x, predicted, group = group, shape = group, linetype = group)) + 
  geom_line(size = .6) + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a cooperative event happens") +
  ggtitle("(a): Cooperative events (Model 6)") + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model 7 (conflictual events, year dummies; no network variables)
ratemodel.7 <- glm(eventDummy ~ 
                     same_annex_factor*yeardummy +
                     log_pop_comp_sender + log_pop_comp_target +
                     log_gdpcap_comp_sender + log_gdpcap_comp_target +
                     log_co2cap_comp_sender + log_co2cap_comp_target + 
                     democracy_comp_sender + democracy_comp_target +
                     nd_gain_sender + nd_gain_target +
                     forest_area_comp_sender + forest_area_comp_target +
                     fossilrents_comp_sender + fossilrents_comp_target +
                     English_sender + English_target +
                     log_trade +
                     log_aid +
                     same_region +
                     same_coalition_noG77 +
                     coalition_member +
                     interv_sender + interv_target +
                     yeardummy,
                   family = binomial(link = logit), data = data_rate_confl)

# Predicted probabilities: 
predict.rm7 <- ggpredict(ratemodel.7, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rm7 <- ggplot(predict.rm7, aes(x, predicted, group = group, shape = group, linetype = group)) + 
  geom_line(size = .6) + 
  geom_point(size=2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a conflictual event happens") +
  ggtitle("Rate Model 7") + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model 8 (conflictual events, year dummies, all variables)
ratemodel.8 <- glm(eventDummy ~ 
                     sender_outdegree + target_indegree +
                     reciprocity + 
                     triad_fof + triad_foe + triad_eoe + triad_eof +
                     sender_similarity + target_similarity +
                     same_annex_factor*yeardummy +
                     log_pop_comp_sender + log_pop_comp_target +
                     log_gdpcap_comp_sender + log_gdpcap_comp_target +
                     log_co2cap_comp_sender + log_co2cap_comp_target +
                     democracy_comp_sender + democracy_comp_target +
                     nd_gain_sender + nd_gain_target +
                     forest_area_comp_sender + forest_area_comp_target +
                     fossilrents_comp_sender + fossilrents_comp_target +
                     English_sender + English_target +
                     log_trade +
                     log_aid +
                     same_region +
                     same_coalition_noG77 +
                     coalition_member +
                     yeardummy,
                   family = binomial(link = logit), data = data_rate_confl)

# Predicted probabilities: 
predict.rm8 <- ggpredict(ratemodel.8, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rm8 <- ggplot(predict.rm8, aes(x, predicted, group = group, shape = group, linetype = group)) + 
  geom_line(size = .6) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a conflictual event happens") +
  ggtitle("(b): Conflictual events (Model 8)") + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


## Figure 5: Effect of group membership on negotiation event frequency over time #####

plots_fig5 <- list(gplot.rm6, gplot.rm8)
margin = theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))
fig_5 <- grid.arrange(grobs = lapply(plots_fig5, "+", margin), nrow = 2, ncol = 1)

ggsave("Figure_5.ps", fig_5, device = "eps", height = 8, width = 6.5, dpi = 300, units = "in")


## Table A9: Results table for appendix #####
htmlreg(list(ratemodel.5, ratemodel.6, ratemodel.7, ratemodel.8), 
        file = "Table_A9.doc", inline.css = FALSE, doctype = TRUE, digits = 3,
        html.tag = TRUE, head.tag = TRUE, body.tag = TRUE, single.row = TRUE)


## Table A10: Robustness checks for the rate model #####

# Model A9 (cooperative events, year polynomial, all variables)
ratemodel.A9 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*year +
                       I(year^2) + I(year^3) +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member,
                     family = binomial(link = logit), data = data_rate_coop)
summary(ratemodel.A9)

# Creating new data on which to run the predictions - the automatic prediction including standard errors does not work
newdata <- expand.grid(year = 1995:2013, same_annex_factor = c("Different annex", "Both non-Annex I", "Both Annex I"))
meanvalues <- sapply(ratemodel.A9$model, mean)
meanvalues <- as.data.frame(t(meanvalues))
meanvalues$same_annex_factor <- NULL
meanvalues$year <- NULL
newdata <- data.frame(newdata, meanvalues)

# Estimating predicted probabilities for the graphs
predict.rmA9 <- predict(ratemodel.A9, type="response", se.fit = FALSE, newdata = newdata)

# Dataset for the graph
gdata.rmA9 <- data.frame(newdata, fit = predict.rmA9)

# Graph
gplot.rmA9 <- ggplot(data=gdata.rmA9, aes(x = year, y = fit, color = same_annex_factor, group = same_annex_factor)) +
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = same_annex_factor), size=.6) + geom_point(aes(shape = same_annex_factor), size=2) + 
  xlab("Year") + ylab("Predicted probability that a cooperative event happens") +
  ggtitle("Rate Model A9") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=7),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A10 (cooperative events, meeting_month (dummies for individual meetings), all variables)
summary(ENBrate_positive$neg_days)
summary(ENBrate_positive$eventTime)
summary(ENBrate_positive$timedummy)
class(ENBrate_positive$meeting_month)

ratemodel.A10 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*meeting_month +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member,
                     family = binomial(link = logit), data = data_rate_coop)
summary(ratemodel.A10)
# takes about 1 hour!

newdata <- expand.grid(meeting_month = data_rate_coop$meeting_month, same_annex_factor = c("Different annex", "Both non-Annex I", "Both Annex I"))
meanvalues <- sapply(ratemodel.A10$model, mean)
meanvalues <- as.data.frame(t(meanvalues))
meanvalues$same_annex_factor <- NULL
meanvalues$meeting_month <- NULL
newdata <- data.frame(newdata, meanvalues)

predict.rmA10 <- predict(ratemodel.A10, type="response", se.fit = TRUE, newdata = newdata)
# takes about 1.5 hours to compute!

gdata.rmA10 <- with(predict.rmA10, 
                      data.frame(newdata, 
                                 fit = fit,
                                 lwr = (fit + 1.96 * se.fit),
                                 upr = (fit - 1.96 * se.fit)))

# Ggplot
gplot.rmA10 <- ggplot(data=gdata.rmA10, aes(x = meeting_month, y = fit, color = same_annex_factor, group = same_annex_factor)) +
  geom_line(aes(linetype = same_annex_factor), size=.6) +
  xlab("Negotiation meeting") + ylab("Predicted probability that a cooperative event happens") +
  ggtitle("Rate Model A10") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text.y = element_text(size=7),
    axis.text.x = element_text(size=5, angle = 90),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A11 (cooperative events, year dummy, no coalitions)
rate_coop_nocoalitions <- subset(data_rate_coop, data_rate_coop$coalition == 0)

ratemodel.A11 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*yeardummy +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       yeardummy,
                     family = binomial(link = logit), data = rate_coop_nocoalitions)
summary(ratemodel.A11)

# Predicted probabilities: 
predict.rmA11 <- ggpredict(ratemodel.A11, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rmA11 <- ggplot(predict.rmA11, aes(x, predicted, group = group, color = group)) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = group), size = .6) + 
  geom_point(aes(shape = group), size = 2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a cooperative event happens") +
  ggtitle("Rate Model A11") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A12 (cooperative events, year dummy, no imputed or completed data)
ratemodel.A12 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*yeardummy +
                       log_population_sender + log_population_target +
                       log_gdpcap_sender + log_gdpcap_target +
                       log_co2cap_sender + log_co2cap_target +
                       democracy_sender + democracy_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_sender + forest_area_target +
                       fossil_rents_sender + fossil_rents_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member +
                       yeardummy,
                     family = binomial(link = logit), data = data_rate_coop)
summary(ratemodel.A12)


# Predicted probabilities: 
predict.rmA12 <- ggpredict(ratemodel.A12, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rmA12 <- ggplot(predict.rmA12, aes(x, predicted, group = group, color = group)) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = group), size = .6) + 
  geom_point(aes(shape = group), size = 2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a cooperative event happens") +
  ggtitle("Rate Model A12") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A13 (conflictual events, year polynomial, all variables)
ratemodel.A13 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*year +
                       I(year^2) + I(year^3) +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member,
                     family = binomial(link = logit), data = data_rate_confl)
summary(ratemodel.A13)

# Creating new data on which to run the predictions - the automatic prediction including standard errors does not work
newdata <- expand.grid(year = 1995:2013, same_annex_factor = c("Different annex", "Both non-Annex I", "Both Annex I"))
meanvalues <- sapply(ratemodel.A13$model, mean)
meanvalues <- as.data.frame(t(meanvalues))
meanvalues$same_annex_factor <- NULL
meanvalues$year <- NULL
newdata <- data.frame(newdata, meanvalues)

# Estimating predicted probabilities for the graphs
predict.rmA13 <- predict(ratemodel.A13, type="response", se.fit = FALSE, newdata = newdata)

# Creating the dataset for the graph
gdata.A13 <- data.frame(newdata, fit = predict.rmA13)

# Ggplot
gplot.rmA13 <- ggplot(data=gdata.A13, aes(x = year, y = fit, color = same_annex_factor, group = same_annex_factor)) +
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = same_annex_factor), size=.6) + geom_point(aes(shape = same_annex_factor), size=2) + 
  xlab("Year") + ylab("Predicted probability that a conflictual event happens") +
  ggtitle("Rate Model A13") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=7),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A14 (conflictual events, meeting_month (dummies for negotiation meetings), all variables)
ratemodel.A14 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*meeting_month +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member,
                     family = binomial(link = logit), data = data_rate_confl)
summary(ratemodel.A14)

# Predicted probabilities: 
predict.rmA14 <- ggpredict(ratemodel.A14, terms = c("meeting_month", "same_annex_factor"))

# Graph:
gplot.rmA14 <- ggplot(predict.rmA14, aes(x, predicted, group = group, color = group)) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = group), size = .6) + 
#  geom_point(aes(shape = group), size = 2) + 
  xlab("Year") + 
  ylab("Predicted probability that a conflictual event happens") +
  ggtitle("Rate Model A14") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text.y = element_text(size=7),
    axis.text.x = element_text(size=5, angle = 90),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A15 (conflictual events, year dummy, no coalitions)
rate_confl_nocoalitions <- subset(data_rate_confl, data_rate_confl$coalition == 0)

ratemodel.A15 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*yeardummy +
                       log_pop_comp_sender + log_pop_comp_target +
                       log_gdpcap_comp_sender + log_gdpcap_comp_target +
                       log_co2cap_comp_sender + log_co2cap_comp_target +
                       democracy_comp_sender + democracy_comp_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_comp_sender + forest_area_comp_target +
                       fossilrents_comp_sender + fossilrents_comp_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       yeardummy,
                     family = binomial(link = logit), data = rate_confl_nocoalitions)
summary(ratemodel.A15)

# Predicted probabilities: 
predict.rmA15 <- ggpredict(ratemodel.A15, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rmA15 <- ggplot(predict.rmA15, aes(x, predicted, group = group, color = group)) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = group), size = .6) + 
  geom_point(aes(shape = group), size = 2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a conflictual event happens") +
  ggtitle("Rate Model A15") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  ylim(0, 0.2) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Model A16 (conflictual events, year dummy, no imputed or completed data)
ratemodel.A16 <- glm(eventDummy ~ 
                       sender_outdegree + target_indegree +
                       reciprocity + 
                       triad_fof + triad_foe + triad_eoe + triad_eof +
                       sender_similarity + target_similarity +
                       same_annex_factor*yeardummy +
                       log_population_sender + log_population_target +
                       log_gdpcap_sender + log_gdpcap_target +
                       log_co2cap_sender + log_co2cap_target +
                       democracy_sender + democracy_target +
                       nd_gain_sender + nd_gain_target +
                       forest_area_sender + forest_area_target +
                       fossil_rents_sender + fossil_rents_target +
                       English_sender + English_target +
                       log_trade +
                       log_aid +
                       same_region +
                       same_coalition_noG77 +
                       coalition_member +
                       yeardummy,
                     family = binomial(link = logit), data = data_rate_confl)
summary(ratemodel.A16)

# Predicted probabilities: 
predict.rmA16 <- ggpredict(ratemodel.A16, terms = c("yeardummy", "same_annex_factor"))

# Graph:
gplot.rmA16 <- ggplot(predict.rmA16, aes(x, predicted, group = group, color = group)) + 
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype = group), size = .6) + 
  geom_point(aes(shape = group), size = 2) + 
  geom_errorbar(aes(ymax=conf.high, ymin=conf.low), width=.02) +
  xlab("Year") + 
  ylab("Predicted probability that a conflictual event happens") +
  ggtitle("Rate Model A16") + 
  guides(color = guide_legend(title = "Group membership")) + 
  guides(linetype = guide_legend(title = "Group membership")) + 
  guides(shape = guide_legend(title = "Group membership")) + 
  scale_x_discrete(breaks=c("1995", "2000", "2005", "2010")) +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))


# Print Table A10
htmlreg(list(ratemodel.A9, ratemodel.A10, ratemodel.A11, ratemodel.A12, ratemodel.A13, ratemodel.A14, ratemodel.A15, ratemodel.A16), 
        file = "Table_A10.doc", inline.css = FALSE, doctype = TRUE, digits = 3,
        html.tag = TRUE, head.tag = TRUE, body.tag = TRUE, single.row = TRUE)


## Figure A13: Plots of further models #####

# Cooperative interactions:
plots_figA13_a <- list(gplot.rmA9, gplot.rmA10, gplot.rmA11, gplot.rmA12)
margin = theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))
fig_A13_a <- grid.arrange(grobs = lapply(plots_figA13_a, "+", margin), ncol = 2)
ggsave("Figure_A13a.pdf", fig_A13_a, device = "pdf", height = 7, width = 15)


# Conflictual interactions:
plots_figA13_b <- list(gplot.rmA13, gplot.rmA14, gplot.rmA15, gplot.rmA16)
fig_A13_b <- grid.arrange(grobs = lapply(plots_figA13_b, "+", margin), ncol = 2)
ggsave("Figure_A13b.pdf", fig_A13_b, device = "pdf", height = 7, width = 15)



###--------------------------------------------------------------------------###
##          Additional test of socialization hypothesis                    #####
###--------------------------------------------------------------------------###

# Use the data_type dataset to create graphs of (mean) reciprocity over time and triads over time, 
# for the whole sample, for each group, and for between-group dyads

library(ggpubr)

## Figure A14: Reciprocity over time #####

graph_recip <- ggplot(data_type, aes(x=year, y=reciprocity, group=same_annex_factor, color=same_annex_factor, fill=same_annex_factor)) + 
  geom_jitter(alpha = 0.8, size = 0.8, shape = 21, color = "white") +
  scale_fill_manual(values = c("palevioletred2", "darkolivegreen3", "orange"),
                    name = "Group membership") +
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4"),
                     name = "Group membership") + 
  stat_summary(aes(x=year, y=reciprocity, color=same_annex_factor), geom="errorbar", fun.data="mean_se",
               size=.3, width=.2) +
  stat_summary(fun = mean, geom = "line") +
  stat_summary(fun = mean, geom = "point", size =.9) + 
  labs(x="Year", y = "Reciprocity") +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))

ggsave(filename="Figure_A14.pdf", graph_recip, width=16, height=8, units="cm")


## Figure A15: Friendship triads over time #####

graph_triad <- ggplot(data_type, aes(x=year, y=triad_fof, group=same_annex_factor, color=same_annex_factor, fill=same_annex_factor)) + 
  geom_jitter(alpha = 0.8, size = 0.8, shape = 21, color = "white") +
  scale_fill_manual(values = c("palevioletred2", "darkolivegreen3", "orange"),
                    name = "Group membership") +
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4"),
                     name = "Group membership") + 
  stat_summary(aes(x=year, y=triad_fof, color=same_annex_factor), geom="errorbar", fun.data="mean_se",
               size=.3, width=.2) +
  stat_summary(fun = mean, geom = "line") +
  stat_summary(fun = mean, geom = "point", size =.9) + 
  labs(x="Year", y = "Triad friend of friends") +
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0.5),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"))

ggsave(filename="Figure_A15.pdf", graph_triad, width=16, height=8, units="cm")



###--------------------------------------------------------------------------###
##          Descriptive statistics                                         #####
###--------------------------------------------------------------------------###

interactions <- data.frame(table(data_type$interaction))
names(interactions)[names(interactions) == "Var1"] <- "Interaction_type"
write.csv(interactions, file = "Table_A1.csv")

senders <- data.frame(table(data_type$sender, data_type$annexi_sender))
senders <- senders[order(-senders$Freq),] 
names(senders)[names(senders) == "Var1"] <- "Country"
names(senders)[names(senders) == "Var2"] <- "AnnexI_member"

senders.annexi <- subset(senders, AnnexI_member == 1)
senders.annexi <- head(senders.annexi, 10)
senders.nonannexi <- subset(senders, AnnexI_member == 0)
senders.nonannexi <- head(senders.nonannexi, 10)

write.csv(senders.annexi, file = "Table_A2a.csv")
write.csv(senders.nonannexi, file = "Table_A2b.csv")

cooperation <- subset(data_type, cooperation == 1)
senders_coop <- data.frame(table(cooperation$sender, cooperation$annexi_sender))
senders_coop <- senders_coop[order(-senders_coop$Freq),] 
names(senders_coop)[names(senders_coop) == "Var1"] <- "Country"
names(senders_coop)[names(senders_coop) == "Var2"] <- "AnnexI_member"

cooperation.annexi <- subset(senders_coop, AnnexI_member == 1)
cooperation.annexi <- head(cooperation.annexi, 10)
cooperation.nonannexi <- subset(senders_coop, AnnexI_member == 0)
cooperation.nonannexi <- head(cooperation.nonannexi, 10)

write.csv(cooperation.annexi, file = "Table_A3a.csv")
write.csv(cooperation.nonannexi, file = "Table_A3b.csv")

conflict <- subset(data_type, cooperation == 0)
senders_confl <- data.frame(table(conflict$sender, conflict$annexi_sender))
senders_confl <- senders_confl[order(-senders_confl$Freq),] 
names(senders_confl)[names(senders_confl) == "Var1"] <- "Country"
names(senders_confl)[names(senders_confl) == "Var2"] <- "AnnexI_member"

conflict.annexi <- subset(senders_confl, AnnexI_member == 1)
conflict.annexi <- head(conflict.annexi, 10)
conflict.nonannexi <- subset(senders_confl, AnnexI_member == 0)
conflict.nonannexi <- head(conflict.nonannexi, 10)

write.csv(conflict.annexi, file = "Table_A4a.csv")
write.csv(conflict.nonannexi, file = "Table_A4b.csv")

senders <- data.frame(table(data_type$sender))
senders <- senders[order(senders$Freq),] 
names(senders)[names(senders) == "Var1"] <- "Country"
View(senders) # data for Table A5

# Parties with zero observations (not in dataset):
# Albania, Andorra, Brunei Darussalam, Guinea, St. Kitts and Nevis, San Marino, Sao Tome and Principe, North Macedonia, Tonga 


# Figure A1: Count of cooperative negotiation interactions over time, by group membership #####
gdata.cooperativeevents <- aggregate(data_type$cooperation, list(Year=data_type$year, Membership=data_type$same_annex_factor), sum)
ENBnew$conflict <- 1 - ENBnew$cooperation
gdata.conflictualevents <- aggregate(ENBnew$conflict, list(Year=ENBnew$year, Membership=ENBnew$annexi.factor3), sum)

names(gdata.cooperativeevents)[names(gdata.cooperativeevents)=="x"] <- "Events"
names(gdata.conflictualevents)[names(gdata.conflictualevents)=="x"] <- "Events"

gdata.cooperativeevents$EventType <- "Cooperative"
gdata.conflictualevents$EventType <- "Conflictual"

gdata.coopconflevents <- rbind(gdata.cooperativeevents, gdata.conflictualevents)
gdata.coopconflevents$EventType <- as.factor(gdata.coopconflevents$EventType)
gdata.coopconflevents$EventType <- factor(gdata.coopconflevents$EventType, levels = c("Cooperative", "Conflictual"))

gplot.coopconfl.events <- ggplot(data=gdata.coopconflevents, aes(x=Year, y=Events, color=Membership, group=Membership)) +
  scale_color_manual(values = c("red", "darkolivegreen4", "sienna4")) + 
  geom_line(aes(linetype=Membership), size=.6) + geom_point(aes(shape=Membership), size=2) + 
  facet_wrap(~EventType, scales = "free_y") + 
  xlab("Year") + ylab("Count of negotiation interactions") +
  ggtitle("Cooperative and conflictual negotiation interactions over time") + 
  theme_bw() +
  theme(
    plot.title = element_text(size=10, face="bold", hjust = 0),
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    legend.text = element_text(size=8),
    legend.title = element_text(size=8, face="bold"),
    legend.position = "bottom")
ggsave(filename="Figure_A1.pdf", plot=gplot.coopconfl.events, width=17, height=8, units="cm")


## Figures A2 - A8: Boxplots of main control variables, by group membership, for years 1995 and 2013 #####

# First exclude the coalitions, and keep only years 1993 and 2013
IVs_1995or2013 <- subset(indep_vars_annual, indep_vars_annual$coalition == 0 & (indep_vars_annual$year == 1995 | indep_vars_annual$year == 2013))

IVs_1995or2013$annexi <- as.factor(IVs_1995or2013$annexi)
levels(IVs_1995or2013$annexi)
levels(IVs_1995or2013$annexi)[levels(IVs_1995or2013$annexi)=="0"] <- "Non-Annex I"
levels(IVs_1995or2013$annexi)[levels(IVs_1995or2013$annexi)=="1"] <- "Annex I"

IVs_1995or2013$yeardummy <- as.factor(IVs_1995or2013$year)

library(dplyr)
library(tibble)
library(ggrepel)

is_maxmin <- function(x) {
  return(x == min(x, na.rm=TRUE) | x == max(x, na.rm=TRUE))
}

is_max <- function(x) {
  return(x == max(x, na.rm=TRUE))
}

is_min <- function(x) {
  return(x == min(x, na.rm=TRUE))
}


dat <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_max(gdpcap_comp), paste0(country, ", ", "\n", round(gdpcap_comp, digits=0), " $"), 
                         ifelse(is_min(gdpcap_comp), paste0(country, ", ", round(gdpcap_comp, digits=0), " $"), as.numeric(NA))))
ggplot(dat, aes(x = annexi, y = gdpcap_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab("GDP PPP per capita (intl. $)") + theme_bw() +
  geom_text(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8))
ggsave(filename="Figure_A2_boxplot_income.pdf", width=16, height=8, units="cm")


dat2 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_maxmin(co2_comp), paste0(country, ", ", "\n", round(co2_comp, digits=0), " kt"), as.numeric(NA)))

ggplot(dat2, aes(x = annexi, y = co2_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab(bquote('Total' ~CO[2]~ 'emissions (kt)')) + theme_bw() +
  geom_text(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8))
ggsave(filename="Figure_A3_boxplot_totalemissions.pdf", width=16, height=8, units="cm")


dat3 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_maxmin(co2cap_comp), paste0(country, ", ", "\n", round(co2cap_comp, digits=2), " t"), as.numeric(NA)))
ggplot(dat3, aes(x = annexi, y = co2cap_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab(bquote(CO[2]~ ' emissions per capita (t)')) + theme_bw() +
  geom_text(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8))
ggsave(filename="Figure_A4_boxplot_emissionscapita.pdf", width=16, height=8, units="cm")


dat4 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(max = ifelse(is_max(democracy_comp), paste0("Several countries, ", "\n", round(democracy_comp, digits=2)), as.numeric(NA)))  %>%
  mutate(min = ifelse(is_min(democracy_comp), paste0(country, ", ", round(democracy_comp, digits=2)), as.numeric(NA)))
ggplot(dat4, aes(x = annexi, y = democracy_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab("Democracy (Freedome House / Polity)") + theme_bw() +
  geom_text(aes(label = max), na.rm = TRUE, hjust = -0.1, size = 2) + 
  geom_text_repel(aes(label = min), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8)) 
ggsave(filename="Figure_A5_boxplot_democracy.pdf", width=16, height=8, units="cm")


dat5 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_maxmin(nd_gain), paste0(country, ", ", round(nd_gain, digits=2), " %"), as.numeric(NA)))
ggplot(dat5, aes(x = annexi, y = nd_gain)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab("ND-Gain index") + theme_bw() +
  geom_text_repel(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8)) 
ggsave(filename="Figure_A6_boxplot_ndgain.pdf", width=16, height=8, units="cm")


dat6 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_maxmin(forest_area_comp), paste0(country, ", ", round(forest_area_comp, digits=2), " %"), as.numeric(NA)))
ggplot(dat6, aes(x = annexi, y = forest_area_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab("Forest area (%)") + theme_bw() +
  geom_text_repel(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8)) 
ggsave(filename="Figure_A7_boxplot_forestarea.pdf", width=16, height=8, units="cm")


dat7 <- IVs_1995or2013 %>%
  group_by(annexi, yeardummy) %>%
  mutate(maxmin = ifelse(is_max(fossilrents_comp), paste0(country, ", ", "\n", round(fossilrents_comp, digits=2), " %"), 
                         ifelse(is_min(fossilrents_comp), paste0("Several countries, ", "\n", round(fossilrents_comp, digits=2), " %"), as.numeric(NA))))
ggplot(dat7, aes(x = annexi, y = fossilrents_comp)) +
  geom_boxplot() + facet_wrap(~year) + xlab("Group membership") + ylab("Fossil rents (% of GDP)") + theme_bw() +
  geom_text(aes(label = maxmin), na.rm = TRUE, hjust = -0.1, size = 2) + 
  theme(
    axis.title.x = element_text(size=8, face="bold"),
    axis.title.y = element_text(size=8, face="bold"),
    axis.text = element_text(size=8)) 
ggsave(filename="Figure_A8_boxplot_fossilrents.pdf", width=16, height=8, units="cm")



